## a wrapper functions for simulating and estimating
## in different ways

## wrapper function for one sim -> estimate
sim_and_est <- function(filename, k,
                        pars, m, T, E, L, phi0) {
  sim <- opmodel(pars, m, T, E, L, phi0)
  dat <- list(T = T, L = L, # E = E[k],
              B_actual = sim$B[,k], R_actual = sim$R[,k],
              C_actual = sim$C[,k], I_actual = sim$I[,k],
              N_actual = sim$N[,k], phi0 = phi0[k],
              p = pars$p[k], w = pars$w[k], M = pars$M[k], q = pars$q[k],
              beta = pars$beta[k], #B0 = B0[k],
              sigR = pars$sigR)
#            lambda = lambda) #sigR = sigR)#, sigI = sigI, sigN = sigN)
  init <- list(# beta = .5,
               B0 = pars$B0[k], #log_sigR = 1+log(sigR),
#               log_sigI = log(pars$sigI), log_sigN = log(pars$sigN),
#               qI = pars$qI, qN = pars$qN,
               eps = rep(0,T+1))
  run_admb(filename, dat, init, re = T, overwrite = T)
}

## Wrapper for handling simulations and running at 
# set parameters first,, then do the sims

## wrapper for simulating right boundary conditions for q
opmodel.depsim <- function(T_dep, dep, reps, pars, m, T, E, L, phi0,
                           errR, errI, errN) {
  q <- rep(0, m)
  for (k in 1:m) {    
    q.opt <- matrix(0, nrow = reps, ncol = 2, dimnames = list(NULL, c("min", "obj")))
    q.init <- 1
    for (i in 1:reps) {
      cat(i); cat(".")
      if (i > 1) q.init <- q.opt[i-1,1]
      q.opt[i,] <- as.numeric(optimize(opmodel.dep, interval = c(.1,100), T_dep = T_dep, dep = dep, k = k,
                                       pars = pars, m = m, T = T, E = E, L = L, phi0 = phi0,
                                       errR = errR, errI = errI, errN = errN))
    }
    q[k] <- q.opt[which(q.opt[,2] == min(q.opt[,2])), 1] # grabs best fit
    cat("\n")
  }
  q
}

## wrapper for simulating the estimation of the real values
sim_estimate <- function(type, N, pars, m, T, E, L, phi0) {
  out <- list()
  
  for (i in 1:N) {
    beta <- rnorm(m, mu, tau)
    q <- rep(1, m)
    pars <- list(beta=beta, mu=mu, tau=tau, B0=B0, p=p, w=w, M=M, sigR=sigR,
                 qI=qI, sigI=sigI, sigN=sigN, qN=qN, q=q)
    dep <- .4 
    reps <- 20
    q <- opmodel.depsim(dep, reps, pars, m, T, E, L, phi0)
    pars$q <- q
    sim <- opmodel(pars, m, T, E, L, phi0)
    out[[i]] <- sim

    if (type == "nlme") {
      filename <- "estimate_nlme"
      dat <- with(sim, list(m = m, T = sim$T, B = t(sim$B), R = t(sim$R), phi0 = sim$phi0))
      inits <- with(sim, list(B0 = sim$B0,
                              log_sigR = rep(log(pars$sigR), m),
                              mu = pars$mu,
                              log_tau = log(pars$tau),
                              eps = rep(0, m)))
      clean_admb(filename)
      write_admb_data(filename, dat)
      write_admb_inits(filename, inits)
      if (i == 1)
        system(paste("admb -r ", filename, sep = ""))
      system(paste(ifelse(.Platform$OS.type == "unix", "./", ""), filename, sep = ""))
      out[[i]] <- c(out[[i]], read_admb_results(filename))
    }
  
    if (type == "bayesian") {
      filename <- "estimate_bay"
      meanR0 <- apply(sim$R * apply(sim$B, 2, function(x) x > median(x)), 2, function(x) mean(x[x>0]))
      dat <- with(sim, list(m = m, T = T, B = t(B), R = t(R), phi0 = phi0, meanR0 = meanR0))
      inits <- with(sim, list(B0 = B0,
                              log_sigR = rep(log(pars$sigR), m),
                              mu = pars$mu,
                              log_tau = log(pars$tau),
                              beta = sim$beta))
      clean_admb(filename)
      write_admb_data(filename, dat)
      write_admb_inits(filename, inits)
      if (i == 1)
        system(paste("admb ", filename, sep = ""))      
      system(paste(ifelse(.Platform$OS.type == "unix", "./", ""), filename, " -mcmc 50000 -mcsave 500 -mcseed 100", sep = ""))
      out[[i]] <- c(out[[i]], read_admb_results(filename))
      mc.out <- read_mcmc_results(filename)
      burnin <- 10
      out[[i]] <- c(out[[i]],
                    list(beta.mc = apply(mc.out[-(1:burnin),(m*2+1):(m*3)], 2, median),
                         B0.mc = apply(mc.out[-(1:burnin),1:m], 2, median),
                         sigR.mc = exp(apply(mc.out[-(1:burnin),(m+1):(2*m)], 2, median)),
                         mu.mc = median(mc.out[-(1:burnin),(3*m+1)]),
                         tau.mc = median(exp(mc.out[-(1:burnin),(3*m+2)]))))
    }
  }
  out
}

## wrapper for doing the simulations w/ both real and bad data
run_full <- function(real = TRUE, est = TRUE, sing = TRUE, run_q = TRUE,
                     filename.sim = "simulate_RE", filename.est = "estimate_nlme",
                     pars, m, T, E, L, phi0) {
  with(pars, {
    if (run_q) {
      dep <- .4 # depletion relative to B0
      reps <- 20
      
      q <- opmodel.depsim(dep, reps, pars, m, T, E, L, phi0)
      pars$q <- q
    }
    
    sim <- opmodel(pars, m, T, E, L, phi0)

    ## fit Bev-Ho individually!

    if (sing) {
      out.sing <- list()
      for (k in 1:m) {
        dat <- with(sim, list(T = sim$T, B = B[,k], R = R[,k],
                              phi0 = phi0[k]))
        inits <- with(sim, list(B0 = B0[k],
                                log_sigR = log(pars$sigR),
                                beta = beta[k])) # sigR should be variable
        out.sing[[k]] <- run_admb(filename.sing, dat, inits, re = FALSE)
      }
    }
    
    ## run meta-analysis method with real data
    if (real) {
      dat <- with(sim, list(m = m, T = sim$T, B = t(sim$B), R = t(sim$R),
                            phi0 = sim$phi0))
      inits <- with(sim, list(B0 = sim$B0,
                              log_sigR = rep(log(pars$sigR), m),
                              mu = pars$mu,
                              log_tau = log(pars$tau),
                              eps = rep(0, m)))
      out.real <- run_admb(filename.est, dat, inits, re = TRUE, overwrite = TRUE)
    }
    
    ## run meta-analysis method with estimated data from random effects model
    if (est) {
      out.sim <- list()
      for (k in 1:m) {
        dat <- list(T = T, L = L, # E = E[k],
                    B_actual = sim$B[,k], R_actual = sim$R[,k],
                    C_actual = sim$C[,k], I_actual = sim$I[,k],
                    N_actual = sim$N[,k], phi0 = phi0[k],
                    p = p[k], w = w[k], M = M[k], q = q[k],
                    beta = beta[k],
                    sigR = sigR)
        ##            lambda = lambda) #sigR = sigR)#, sigI = sigI, sigN = sigN)
        
        init <- list(# beta = .5,
                     B0 = rnorm(1, 40, 2),
                     ##             log_sigR = rnorm(1, log(sigR), .4),
                     ##             log_sigI = rnorm(1, log(sigI), .4),
                     ##             log_sigN = rnorm(1, log(sigN), .4),
                     u_I = 1.0,
                     u_N = 1.0,
                     qI = qI,
                     qN = qN,
                     eps = rep(0, T+1))
        
        out.sim[[k]] <- run_admb(filename.sim, dat, init, re = T, overwrite = TRUE)
      }
      
      out.BR <- lapply(out.sim, function(x) data.frame(B = x$B, R = x$R))
      out.BR <- data.frame(.id = rep(1:m, each = T+1), do.call(rbind, out.BR))
      
      out.BR[,1] <- as.numeric(out.BR[,1])
      out.B <- unstack(out.BR, terms(B ~ .id, keep.order = TRUE))
      out.R <- unstack(out.BR, terms(R ~ .id, keep.order = TRUE))
      
      dat <- with(sim, list(m = m, T = sim$T, B = t(out.B), R = t(out.R),
                            phi0 = sim$phi0))
      inits <- with(sim, list(B0 = sim$B0,
                              log_sigR = rep(log(pars$sigR), m),
                              mu = pars$mu,
                              log_tau = log(pars$tau),
                              eps = rep(0, m)))
      ## sigR should be variable
      
      out.est <- run_admb(filename.est, dat, inits, re = TRUE, overwrite = TRUE)
    }
    
    ## grab results from
    result <- data.frame(beta = beta,
                         beta_real = with(out.real, eps*exp(log_tau)+mu),
                         beta_est = with(out.est, eps*exp(log_tau)+mu),
                         beta_sing = unlist(lapply(out.sing, function(x) x$beta)),
                         mu_real = out.real$mu, mu_est = out.est$mu,
                         tau_real = exp(out.real$log_tau), tau_est = exp(out.est$log_tau),
                         posdef = unlist(lapply(out.sim, function(x) x$posdef)))
    return(result)
  })
}

## run operating model with fast q optimization
run_opmodel <- function(pars, m, T, E, L, phi0) {
  sim <- opmodel(pars, m, T, E, L, phi0)
  errR <- sim$errR; errI <- sim$errI; errN <- sim$errN

  q.temp <- rep(0, length = m)
  for (k in 1:m) {
    q.temp[k] <- optimize(opmodel.dep, interval = c(.1,50), T_dep = T_dep, dep = dep, k = k,
                          pars = pars, m = m, T = T, E = E, L = L, phi0 = phi0,
                          errR = errR, errI = errI, errN = errN)$minimum
  }
  
  pars$q <- q.temp
  opmodel(pars, m, T, E, L, phi0, errR, errI, errN)
}

run_estmodel <- function(sim, pars) {
  est <- list()
  for (k in 1:m) {
    clean_admb("simulate_RE")
    
    dat <- list(T = T, L = L, 
                B_actual = sim$B[,k], R_actual = sim$R[,k],
                C_actual = sim$C[,k], I_actual = sim$I[,k],
                N_actual = sim$N[,k], phi0 = phi0[k],
                p = p[k], w = w[k], M = M[k], q = q[k],
                beta = pars$beta[k])
    ##              sigR = sigR,
    ##              sigN = sigN)
    
    inits <- list(# beta = beta[k],
                  B0 = B0[k],
                  log_sigR = log(sigR),
                  eps = rep(0, T))
    ##                log_sigI = log(sigI),
    ##                log_sigN = log(sigN),
    ##                qI = qI,
    ##                qN = qN,
    ##                eps = rep(0, T)
    
    write_admb_data("simulate_RE", dat)
    write_admb_inits("simulate_RE", inits)
    
    #system("admb -r simulate_RE")
    system("./simulate_RE")
    cat(paste("\n", k, "\n"))
    
    est[[k]] <- read_admb_results("simulate_RE")
  }
  est
}

run_meta <- function(est, sim, pars, method = c("nlme", "bay", "prof"), real = FALSE) {
  if (!real) {
    est.BR <- lapply(est, function(x) data.frame(B = x$B, R = x$R))
    est.BR <- data.frame(.id = rep(1:m, each = T+1), do.call(rbind, est.BR))
    
    est.BR[,1] <- as.numeric(est.BR[,1])
    est.B <- unstack(est.BR, terms(B ~ .id, keep.order = TRUE))
    est.R <- unstack(est.BR, terms(R ~ .id, keep.order = TRUE))
  }
  if (real) {
    est.B <- sim$B
    est.R <- sim$R
  }
  
  if (method[1] == "nlme") {
    clean_admb("estimate_nlme")
    dat <- with(sim, list(m = m, T = T, L = L, B = t(est.B), R = t(est.R),
                          phi0 = phi0))
    inits <- with(sim, list(B0 = B0,
                            log_sigR = rep(log(pars$sigR), m),
                            mu = pars$mu,
                            log_tau = log(pars$tau),
                            eps = rep(0, m)))

    write_admb_data("estimate_nlme", dat)
    write_admb_inits("estimate_nlme", inits)
    
    #system("admb -r estimate_nlme")
    if (.Platform$OS.type == "unix")
      system("./estimate_nlme")
    if (.Platform$OS.type == "windows")
      shell("estimate_nlme")

    out <- read_admb_results("estimate_nlme")
  }
  
  if (method[1] == "bay") {
    clean_admb("estimate_bay")
    dat <- with(sim, list(m = m, T = T, L = L, B = t(est.B), R = t(est.R),
                          phi0 = phi0, meanR0 = meanR0))
    inits <- with(sim, list(B0 = B0,
                            log_sigR = rep(log(pars$sigR), m),
                            beta = pars$beta,
                            mu = pars$mu,
                            log_tau = log(pars$tau)))
    
    write_admb_data("estimate_bay", dat)
    write_admb_inits("estimate_bay", inits)

    #system("admb estimate_bay")
    system("./estimate_bay -mcmc 50000 -mcsave 200")

    out <- read_admb_results("estimate_bay")
    out$mcmc <- read_mcmc_results("estimate_bay", inits)

    out$B0_mc <- out$mcmc[,grep("B0", names(out$mcmc))]
    out$beta_mc <- out$mcmc[,grep("beta", names(out$mcmc))]
    out$log_sigR_mc <- out$mcmc[,grep("log_sigR", names(out$mcmc))]
    out$mu_mc <- out$mcmc[,grep("mu", names(out$mcmc))]
    out$log_tau_mc <- out$mcmc[,grep("log_tau", names(out$mcmc))]
  }
  
  if (method[1] == "prof") {
    
  }

  out
}
